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We have performed large-scale Monte Carlo simulations on a two-dimensional generalized Ashkin- 
Teller model to calculate the thermodynamic properties in the critical region near its transitions. 
The Ashkin- Teller model has a pair of Ising spins at each site which interact with neighboring spins 
through pair-wise and 4-spin interactions. The model represents the interactions between orbital 
current loops in Cu02-plaquettes of high-T c cuprates, which order with a staggered magnetization 
M s inside each unit-cell in the underdoped region of the phase diagram below a temperature T* (x) 
which depends on doping. The pair of Ising spins per unit-cell represent the directions of the currents 
in the links of the current loops. The generalizations are the inclusion of anisotropy in the pair- 
wise nearest neighbor current-current couplings consistent with the symmetries of a square lattice 
and the next nearest neighbor pair- wise couplings. We use the Binder cumulant to estimate the 
correlation length exponent v and the order parameter exponent f3. Our principal results are that in 
a range of parameters, the Ashkin- Teller model as well as its generalization has an order parameter 
susceptibility which diverges as T — > T* and an order parameter below T* . Importantly, however, 
there is no divergence in the specific heat. This puts the properties of the model in accord with the 
experimental results in the underdoped cuprates. We also calculate the magnitude of the "bump" 
in the specific heat in the critical region to put limits on its observability. Finally, we show that the 
staggered magnetization couples to the uniform magnetization Mq such that the latter has a weak 
singularity at T* and also displays a wide critical region, also in accord with recent experiments. 

PACS numbers: 74.20.Rp, 74.50.+r, 74.20.-z 
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I. INTRODUCTION 

It has been proposed^* 2 - that the properties of the 
cuprate compounds are controlled by the onset of a time- 
reversal and inversion violating order parameter below a 
temperature T = T*{x), which depends on the doping 
x. T*(x) — * for x — > x c in the superconducting range 
of compositions, thus defining a quantum critical point. 
The quantum critical fluctuations associated with the 
breakup of the specific order proposed have been shown 3 
to be of the scale-invariant form hypothesized to lead to 
a Marginal Fermi Liquid^, which explains the anomalous 
transport properties of these compounds. T*{x) is identi- 
fied with the observed onset of the pseudogap properties 
in the cuprates. 

A major difficulty in accepting these ideas is that there 
is no observed specific heat divergence near T* (x) in any 
cuprate. On the other hand, there now exists signifi- 
cant evidence for long-range order with a spatial symme- 
try consistent with orbital currents of the form shown in 
Fig. [1] in three different families of cuprates^liS which 
have been investigated so far. There is also evidence 
of a weak singularity at T*(x) in the uniform magnetic 
susceptibility. 

In view of this situation, it is important to investi- 
gate whether or not the proposed models for these novel 
broken symmetries are consistent simultaneously with 
long-range order without an observable signal in the spe- 
cific heat in the measurements made hitherto, and also 
whether it does give rise to observable features in the 



uniform magnetization induced by an external magnetic 
fields 

The particular form of proposed hidden order is one of 
spontaneously generated fluxes in the O-Cu-O plaquettes 
of the CmC?2 unit cell such that currents flow in two oppo- 
sitely directed loops in each unit-cell, as depicted for one 
of the four possible domains in Fig. Q] (See also Fig. 1 
of Refs. I^fioh. The staggered orbital magnetic moments 
within each CuOi unit cell repeats from unit cell to unit 
cell so that the translational symmetry of the lattice re- 
mains unaltered. These circulating current patterns are 
generated by a nearest-neighbor repulsion V between Cu 
and O-atoms in the CwC^-sheets. The effect of such a re- 
pulsive I^-term has been extensively investigated in ID 
CuO-chains, where it has been shown to drive charge- 
transfer instabilities and superconductivit y 11 : 12 ' . In 
Ref. Ihll the existence of current-loop ordering was not 
confirmed, but the work was carried out on a truncated 
effective t— J model of 8 Cu sites. Moreover, the ground 
state was spin-polarised with finite momentum, which 
would not be representative of the large-scale physics of 
interest in the system. The truncation of the Hilbert 
space used in Ref. lij furthermore requires so large val- 
ues of onsite Coulomb repulsion on oxygen sites that it 
is probably outside the parameter regime of the high-T c 
cuprates 15 . This motivated the authors of Ref. [l^to un- 
dertake a large-scale study of the issue of current-loop 
ordering on much larger systems using the full three- 
band model of the CuOi planes via variational Monte 
Carlo simulations. These authors find clear evidence for 
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current-loop ordering. Other types of current-patterns 
and charge-fluctuations are also possible^ ' 17 ' 18 . 




O sites fl| Cu sites 



FIG. 1: (Color online) The circulating current phase 0/j— . 
The Cu sites are red circles, O sites are blue. The unit cell is 
shown by the dashed square. A staggered magnetic moment 
pattern within each unit cell that repeats from unit cell to 
unit cell (the curl of the directed circles) is indicated. The 
currents J x and J y represent the horizontal end vertical cur- 
rents, respectively, to be used in the derived effective model, 
Eq. [T] below. Physically, they represent the coherent parts of 
the orbital fermionic currents in the problem. 



II. FUNDAMENTALS 

In this section, we present the effective model of fluctu- 
ating orbital currents we study in this paper, along with 
the definitions of the thermodynamic quantities we com- 
pute, as well as some remarks on the critical exponents 
of the problem with emphasis on the particular status of 
the specific heat exponent of the problem at hand. 



A. Model 

The effective model we perform Monte Carlo simula- 
tions on, has been derived from a microscopic description 
of the CuCVplanes of high-T c cuprates elsewher o 10 i 19 . It 
turns out to be a generalization of the model initially 
proposed to describe the statistical mechanics of loop 
current order— ^ in which some terms allowed by sym- 
metry were omitted. The action S is written on the form 
S = Sc + Sq, where Sc is the classical piece of the ac- 
tion, and Sq is part of the action that is needed in the 
quantum domain of the theory. In this paper, we will fo- 
cus on discussing the effects of thermal fluctuations, and 
we will therefore not need .Sq. The classical part of the 



action, Sc, is given b y 10 i 19 

S C = -0 £ i K ' JrJr'+Ky J V r J V v ,) 
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Here, (r, r') and ((r, r')) denote nearest-neighbor and 
next-nearest-neighbor summations, respectively. = 
1/T where T is temperature, and we work in units where 
Boltzmann's constant ks — 1. We will only consider the 
directions ± of the current variables J x ' v , assuming as 
in other similar two-dimensional models that their am- 
plitudes are smoothly varying with temperature and do 
not determine the critical properties. Note that there is 
always also a current in the O — O links whose magnitude 
is equal to that of J x which has the same magnitude as 
J y . Therefore, no current flows out of any O — Cu — O tri- 
angular plaquette. Due to the restriction that no current 
flows out of any O — Cu — O plaquette, there is no need 
to specify the O — O currents. The variables J x , J y are 
then the same as the a = ±1 and r = ±1 Ising variables 
introduced earlier 3 . Fluctuations (J* — > — J x ,jv — » J$) 
corresponds to going from the depicted current pattern 
(Fig. rjj to a new one which is obtained by a counterclock- 
wise rotation by ir/2, (J x — ► J x ,Jj! — » — J]!) corresponds 
to clockwise rotation of n/2, and (J* — ► — J x , J y — ► — ) 
to a rotation of n. 

If one ignores the next-nearest neighbor terms and 
takes K x = K y , one gets the Ashkin- Teller (AT) models, 
for which several exact results are known^i asymptoti- 
cally close to the phase transition lines. However, since 
the currents are bond-variables, one necessarily has an 
anisotropy in the nearest neighbor interaction o 10 ' 19 , such 
that for r — r' = ±x, K x = Ki and K y = K t , whereas 
when r — r' = ±y, K x — K t and K y = Ki. It is im- 
portant to investigate whether this anisotropy is an ir- 
relevant perturbation. We will in the following denote 
the anisotropy by the parameter A = K t /Ki. Sim- 
ilarly, it is interesting to investigate the effect of the 
next-nearest neighbor interaction given by the parameter 
K xv = K xy when r _ r ' = ±(x + y ) an d = ~K x v 

when r — r' = ±(x — y). 

Let us comment briefly on the terms appearing to quar- 
tic order, most of which cither arc constants or renormal- 
ize the quadratic piece of the action. Note that four Ising 
variables of two distinct species all located on one single 
lattice site, simply contribute a constant to the action. 
If we now limit ourselves to terms that have four J-fields 
distributed on two nearest-neighbor lattice sites, only two 
distinct possibilities exist. Firstly, we may have a term 
with three J's on one lattice site and one J on a nearest- 
neighbor site. This merely represents a renormalization 
of the quadratic couplings. Secondly, we may have two 
J's on one lattice site and another two on a nearest neigh- 
bor lattice site. Unless there are two distinct species of 
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J's on each of the lattice sites, such a term will represent 
a constant contribution to the action. If the J's on each 
lattice site are of distinct species, the term will be of the 
AT-form, as written above. We will ignore terms that 
have J-fields distributed on three or four distinct lattice 
sites, such as for instance plaquette terms, as these are 
generated by much higher order term o 10 ' 19 . 



B. Thermodynamic quantities 

In this paper, we calculate the evolution of the specific 
heat, the staggered orbital magnetic moment as well as 
the susceptibility of the staggered orbital magnetic mo- 
ment as we vary K4 in Eq. Q] We also perform finite-size 
scaling on the magnetization and the Binder cumulant 
(see below). The specific heat C v is given by 



a = -l((5 c -(5 c )) 2 ). 



(2) 



Considering Fig. [21 we see that we may define a pseudo- 
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FIG. 2: (Color online) An illustration of the pseudo- "spin" 
S = ( Jr , Jr ) we use to compute the staggered order parame- 
ter and its susceptibility, Eqs. [3] and [4] 

"spin" S on each lattice given by S r = {J*,Jj>). The 
various states of the system are then described by a 4- 
state clock pseudospin S r = (±1, ±1) on a 2-dimcnsional 
square lattice. We define the staggered order parameter 
in the standard way it would be defined for a clock model, 
namely 
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where m a = ^ r ^r*' a ^ ( x > f )• The susceptibility of this 
staggered order parameter is given by 
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((m x ) 2 + {m y f) - (y/(m x ) 2 + {mv) 2 ) 2 .(4) 



We will contrast the singularities in these quantities with 
the evolution of the anomaly in the specific heat as the 
parameter K4 is varied. While the above staggered mo- 
ment does not couple linearly to an external uniform 
magnetic field, it couples to a field-induced uniform mag- 
netic moment via a quartic term in the free energy. The 
field-induced uniform magnetization must therefore have 
a non-analytic behavior across the phase transition where 
the staggered magnetization associated with the ordering 
of the orbital currents sets in. We will return to this point 
in Section [Tvl 

For the purposes of extracting the critical exponent v, 
we consider the Binder cumulant, defined by 
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where m 2 — (m x ) 2 + (m v ) 2 , corresponding to the magne- 
tization order parameter (|m|), whose critical exponent 
P is given in Eq. [6] for the AT-model 2 ^. In the ordered 
phase, G = 1. For an A-component order parameter, 
G = (N + 2)/N in the disordered phase. In our case, 
therefore, G will exhibit a rise from 1 to 2 as the sys- 
tems disorders. When computing this quantity for dif- 
ferent L and plotting it as a function of T, the curves 
should in principle cross at the same point, thus defin- 
ing T c . On the other hand, plotting it as a function of 
L 1 /' J \(T - T c )/T c )\, all the curves will collapse on top of 
each other. By adjusting v to get data-collapse, one ob- 
tains the correlation length exponent. Furthermore, the 
order-parameter exponent /? is obtained from the mag- 
netization M s for various system sizes by considering the 
quantity L^l v ' M s and adjusting (3 and v so as to obtain 
data-collapse when plotting this quantity as a function 
of L l / V \(T - T c )/T c )\. 



C. Critical exponents 

Note that although the K x and K y couplings between 
the two different types of Ising fields in this model are 
anisotropic ! 10 ! 19 ' 22 , there is only one (doubly degenerate 
Ising) phase transition in the system for K xy — 0; K4 = 
0. Hence, as the four-spin coupling K4 is changed from 
0, the Ising critical point evolves into a single phase- 
transition line with non-universal critical exponents 21 . 
In particular, the specific heat exponent a becomes nega- 
tive, with the transition line itself being a selfdual critical 
line^i. In this sense, the model is similar to an isotropic 
AT model, where the exact result for the critical expo- 
nents are known, and given by2i 



2- 2?/ 1 
a = : a = — , 

3- 2y' H 8 \i-2y 



(G) 



From this, we deduce the susceptibility exponent 7 = 
14/3 and the correlation length exponent v = 8/3 from 
standard scaling relations. Note that the ratios j/v = 
7/4 and f3/v =1/8 are universal and independent of y. 
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(It is also interesting to note that the anomalous scaling 
dimension r\ = 1/4 and the magnetic field exponent 5 = 
15, precisely as in the 2D Ising model). Here y = 2/j/it 
and cos(jti) = [e 4Ki / T - - l]/2 21 . Hence, for K 4 < 0, we 
have 7r/2 < \i < 2tt/3, such that 1 < y < 4/3. 

These exponents are plotted in Fig. [31 The most ex- 
treme deviation from the 2D Ising values a = 0, = 
1/8, 7 = 7/4, v = 1 is given by the case K4 — > — oo, y = 
4/3, where a = -2, /3 = 1/4, 7 = 7/2, v = 2. Note the 
increase of 7 and ^, (which implies a weak increase in 
f3 for increasing — K4, due to the proportionality factors 
14 and 8 given below Eq. [6]), while we have a substan- 
tial reduction of a to negative values as —K4 increases. 
This is traceable to the numerator 2 — 2y in a compared 
to the numerator 2 — y in /3, 7, and ^ (while 77 and 5 
are independent of y). Hence, the specific heat exponent 
stands out as very special in the model Eq. [1] This fact 
is by far the single most dramatic difference between the 
critical behavior of Eq. Q] and the 2D Ising model. The 
-K4-term with K4 < simultaneously suppresses singu- 
larities in the specific heat, and enhances singularities 
both in the susceptibility corresponding to the staggered 
orbital magnetization of Fig. [T]and in the one associated 
with a field-induced uniform magnetization (see Section 

IE}. 
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FIG. 3: (Color online) Critical exponents a, j3, and 7 from the 
Ashkin- Teller model, as a function of the four-spin coupling 
PcKi < O^i. In this parameter range, we have —2 < a < 0, 
1/8 < P < 1/4, 7/4 < 7 < 7/2, and 1 < v < 2. 



III. MONTE CARLO RESULTS 

The Monte Carlo computations were performed us- 
ing the standard single-spin update Metropolis-Hastings 
algorith m 23 i 24 , making local updates of the Ising-fields 
J x and J y , as well as local updates of the composite 
Ising-field J x J y at each lattice site. The system-grid is 
defined by two 2-dimensional subgrids, one for each Ising- 
field, and the local updates were performed for all points 
on the grid. All the Ising-fields on both subgrids were 



initially set to 1. We started all simulations at the high- 
temperature end, and discarded the first 100000 sweeps 
for the purpose of initially thermalizing the system. Af- 
ter that, measurements were made for every 100 sweeps. 
The system sizes that were considered were L x L with 
L = 64, 128, 256, 512. For each value of T, we ran up to 
3 ■ 10 6 MC sweeps for L = 64, 128, 256 and sampled the 
system for every 100 MC sweeps over the lattice, while we 
used 5 • 10 6 MC sweeps for L = 512 and sampled the sys- 
tem for every 150 MC sweeps over the lattice. We have 
checked that satisfactory convergence is well established 
by the time we get to system sizes of L = 512, and we 
therefore largely present results for these largest systems 
only, apart from Fig. [5] and the finite-size scaling re- 
sults that will be presented for the Binder-cumulant(see 
below). In all simulations, we have set K\ — 1.0, such 
that all other couplings are measured relative to this pa- 
rameter. In these units, the critical temperature T c of 
the system for A = 1.0, K xy = Q,K^ = is given by 
T c = 2/ln(l + V2) fa 2.27. This sets the scale of the 
critical temperatures in the plots we will show below. 



A. Specific heat 

Let us first investigate what effect K xy has on the log- 
arithmic singularity of the 2D Ising model. In Fig. [U 
we show the specific heat for A = 1.0 and K4 = 0, 
upon varying K xy = 0.0,0.1,0.2,0.3. We have limited 
the variations in K xy because it can be shown in mean- 
field calculations that the order parameter changes the 
translational symmetry for large enough K xy and a di- 
agonal "striped" order is favored. It is seen that the 
K xy term in this parameter range leaves the logarithmic 
singularity of the anisotropic double-Ising model (Eq. Q] 
with K xy = 0, K4 = 0) unaltered, only the amplitude of 
the singularity is changed. 

We now investigate the effect of four-spin interactions 
oc K4 . We will only consider negative values of K4 in this 
paper. Then the four-spin term tends to promote a non- 
uniform ground state with antiferromagnetic ordering in 
the composite variable J x J y , thus frustrating the Ising 
terms in Eq. [TJ It is known from the phase diagram of 
the AT model 2 ^ that the ordered phase has a different 
symmetry in the regions — 1 < K4/K1 < l,Ki/Ki < 

— 1 and K4/K1 > 1. The region of special interest is 

— 1 < K4/K1 < in which the AT model has a self-dual 
line of critical point o 25 ' 26 . This is consistent with the 
microscopic model, which may exhibit a negative sign of 
the four-spin interaction term. 

We first consider the case of isotropic Ising coupling 
Ki = K t , i.e. A = 1.0, next-nearest neighbor coupling 
K xy = 0.0, and increasing \K^\. We use this case for 
reference, as this parameter set represents the standard 
isotropic AT modet^ii. The results for the specific heat 
are shown in Fig. \5\ The logarithmic specific heat of the 
Ising model disappears to be replaced by a bump whose 
extent in T increases as \K^\ increases. This is consistent 
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FIG. 4: (Color online) Specific heat as a function of temper- 
ature T for the classical part of the model in Eq. [1] with A = 
1.0 and K 4 = 0.0, for various values of K xy = 0.0, 0.1, 0.2, 0.3, 
and system size L = 512. The amplitude of the logarithmic 
specific heat of the Ising model (K xy = 0), is enhanced as 
K xy increases, but the anomaly remains logarithmic. The 
critical temperature of the 2D pure Ising model is given by 
T c = 2/ln(l + v2) ~ 2.27 in units where Boltzmanns con- 
stant fcs = 1. Note also that for this set of parameters, K xy 
hardly alters T c of the model with K xy = 0. 



with the asymptotic critical exponents 
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FIG. 5: (Color online) Specific heat as a function of temper- 
ature T for the classical part of the generalized AT model 
Eq. [TJ with A = 1.0 and K xy — 0.0, for various values of 
K 4 = 0.0, -0.1, -0.25, -0.5, and system size L = 512. The 
vertical scale is in units of fcs /unit-cell. The logarithmic spe- 
cific heat singularity of the Ising model (K 4 = 0), is elimi- 
nated and replaced by a bump whose width increases as \Ka\ 
increases. The arrow in the lower right panel indicates T c as 
obtained from the peak in the susceptibility \ s - 

In Fig. [SI we investigate how well these results are 
converged when increasing the system size through the 
values L = 64, 128, 256, 512. It is seen that the results 
appear well converged when L has reached 256, in par- 
ticular the double-peak structure in Cy that is present 
for small system sizes disappears upon increasing L. In 



contrast to the Binder-cumulant (see below), we have not 
attempted a data collapse of the specific heat by trying 
a scaling form C V (T, L) = L a l v C± (L 1 ^ (T — T c )/T c ) and 
adjusting a to obtain data-collapse. The reason is that 
we anticipate a negative specific heat exponent, such that 
corrections to the above scaling form will be large, thus 
preventing data collapse. Even for positive a, it is well- 
known that corrections to scaling are substantial for the 
specific heat. This simply means that the specific heat 
by itself oddly enough is not a very useful quantity from 
which to extract precise values of a in Monte-Carlo com- 
putations on practical system sizes. Other techniques 
are required for this, see e.g. Ref. : 2f. However, the 
main point of the present paper is not to determine a 
precise value of a numerically, but rather to demonstrate 
(including all corrections to scaling) that a striking sup- 
pression of the prominent logarithmic singularity of the 
2D Ising model takes place as \K 4 \ is increased. Fig. [5] 
clearly shows that the suppression is not a finite-size ar- 
tifact. Note in particular that the relative height of the 
bump in Cy for non-zero \K^\ is suppressed compared to 
the Ising-singularity as L increases. 
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FIG. 6: (Color online) Specific heat as a function of tempera- 
ture T for the classical part of the generalized AT model Eq. 
[TJ with A = 1.0 and K xy = 0.0, for various values of K 4 = 
0.0, -0.1, -0.25, -0.5, and system size L = 64, 128, 256, 512. 
The vertical scale is in units of ks /unit-cell. The logarithmic 
specific heat singularity of the Ising model (if 4 = 0), is elimi- 
nated and replaced by a bump whose width increases as \Ka\ 
increases. Note how the double-bump in Cv , which is present 
at smaller system sizes, disappears when L is increased. When 
L — 512, the results appear to be well converged. 

We next consider the effect of increasing the anisotropy 
(A < 1), such as to weaken the ordering in each of the 
J y (r)- and J x (r) Ising fields. Note, however, that because 
the anisotropy introduced is equal for both of the Ising 
fields (only the direction of the anisotropy is changed) 
the model only has one single critical point even in the 
absence of a ^-coupling. The model is then merely 
two copies of one and the same anisotropic 2D Ising 
model. However, an increase in \K± \ is expected to have a 
stronger effect for A < 1.0 than when A = 1.0 due to the 
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weaker ordering and reduced critical temperature. The 
bump in the specific heat is then accordingly smoother 
as seen in Fig. [7] compared to Fig. [5] 




T 



FIG. 7: (Color online) Specific heat as a function of temper- 
ature T for the classical part of the generalized AT model 
Eq. [TJ with A = 0.5 and K xy = 0.0, for various values of 
K4, = 0.0, —0.1, —0.25, and system size L — 512. The vertical 
scale is in units of fes /unit-cell. Compared to the case shown 
in Fig. with A = 1.0, precisely the same trends are seen in 
the evolution of the anomaly as the AT coupling \K±\ is in- 
creased, only slightly more pronounced. The arrow indicates 
T c as obtained from the peak in the susceptibility Xs- 

We now repeat the above computations for A = 1.0 
with K xy = 0.1, 0.2 and 0.3. This coupling tends to frus- 
trate the Ising ordering, since a large K xy tends to pro- 
mote striped order due to the diagonal anisotropy (repre- 
sented by a change of sign in K xy upon 7r/2 rotations of 
next-nearest neighbor vectors). It is of interest to see how 
the presence of K xy affects the introduction of the AT 
coupling K4, Naively, since the coupling K xy promotes 
striped order and frustrates the uniform order promoted 
by K x , K y , we would expect that the suppressed anoma- 
lies are pushed to lower temperatures as K xy is increased. 
In Figs. [51 [5J and ITTJ1 we show the specific heat for the 
same sets of parameters as in Fig. [5j except that now 
K xy = 0.1, 0.2, 0.3, respectively. 

We see that the effect of K xy is to increase the sharp- 
ness of the bump in the specific heat, while the effect of 
K4 again is to widen the bump (in the presence of K xy ). 
We also see that the anomalies that remain are pushed 
slightly downwards in temperature compared to the case 
K xy — 0, cf. the results of Fig. \5[ The change is how- 
ever only minor for the cases K xy = 0.1 and K xy = 0.2, 
consistent with the weak suppression of the critical tem- 
perature we found upon increasing K xy at K4 = in Fig. 
[4j The conclusion we draw from these computations is 
that the singularity of the specific heat of the Ising case 
is removed by the coupling K4 is included. The resulting 
bump in the specific heat becomes sharper for increasing 
K xy at finite K 4 . 

Finally, we consider the most general case of 
anisotropic Ising coupling A = 0.5 and finite K xy = 0.3, 



K\ = 0.06 ^ 
Ki = -0.10 
7 - Ki = -0.25 » 
Ki = -0.50 □ 
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FIG. 8: (Color online) Specific heat as a function of temper- 
ature T for the classical part of the generalized AT model 
Eq. m, with A = 1.0 and K xy = 0.1, for various values of 
K 4 = 0.0, -0.1, -0.25, -0.5, and system size L = 512. The 
arrow indicates T c as obtained from the peak in the suscepti- 
bility Xs ■ The vertical scale is in units of ka /unit-cell. 
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FIG. 9: (Color online) Specific heat as a function of temper- 
ature T for the classical part of the generalized AT model 
Eq. [1] with A = 1.0 and K xy = 0.2, for various values of 
K 4 = 0.0, -0.1, -0.25, -0.5, and system size L = 512. The 
arrow indicates T c as obtained from the peak in the suscepti- 
bility Xs- The vertical scale is in units of ks/ unit- cell. 



as I if 4 1 is increased, shown in Fig. [TTj It is clear from Fig. 
[IT] that the introduction of anisotropy A = K t /Ki = 0.5 
widens the width of the bump in the specific heat. This 
is easily understood, since increasing anisotropy implies 
that the magnitude of K4 relative to the Ising couplings 
in the problem will increase. The effect of a given in- 
crease in K4 is therefore more strongly felt. Moreover, 
as in the isotropic case, the anomalies are pushed down 
in temperature compared to the case K xy — 0, cf. the 
results of Fig. [7] 

Concluding this section on the results for the specific 
heat, we mention that we have also, at the early stages 
of this work, performed a rather rudimentary compara- 
tive study of the specific heat anomaly in the 2D Ashkin- 
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FIG. 10: (Color online) Specific heat anomaly as a function 
of temperature T for the classical part of the generalized AT 
model Eqs. {TJ, with A = 1.0 and K xy = 0.3, for various 
values of K 4 — 0.0, —0.1, —0.25, —0.5, and system size L = 
512. The vertical scale is in units of fcs /unit-cell. 
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FIG. 11: (Color online) Specific heat anomaly as a function 
of temperature T for the classical part of the generalized AT 
model Eq. [T] with A = 0.5 and K xy = 0.3, for various values 
of K 4 = 0.0,-0.1,-0.25, and system size L = 512. The 
vertical scale is in units of fcs /unit-cell. 



Teller model and the 2DXY continuous rotor model with 
a 4-fold symmetry breaking term, on lattice sites up to 
L = 32. This numerics is insufficient to draw any conclu- 
sions about the fluctuation spectrum on the disordered 
side of the transition, close to the transition, as the sym- 
metry breaking field becomes small. That is, the simula- 
tions per se do not allow us to conclude anything about 
the perturbative relevance or irrelevance of the symme- 
try breaking term. What we have been able to confirm, 
is that the specific heat anomaly of the 2D Ashkin- Teller 
model is indistinguishable from the 2DXY continuous 
rotor model with a symmetry breaking term, provided 
the symmetry breaking term is large. 



B. M s , Xs, and the critical exponents v and (3 

Let us now study the order parameter and suscepti- 
bility of the order parameter, M s and Xs, Eqs. [3] and [J] 
We have first chosen parameters A = 1.0, K xy — 0, and 
varied K4, for which the evolution of the specific heat 
anomaly is shown in Fig. [5] The results for M s and 
Xs are shown in Figs. [T2] and [131 respectively. We see 
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FIG. 12: (Color online) The staggered order parameter, Eq. 
131 as a function of temperature T for the classical part of the 
generalized AT model Eq. Q] with A = 1.0 and K xy = 0.0, 
for various values of K 4 = 0.0, —0.1, —0.25, —0.5, and system 
size L = 512. 

that the staggered magnetization retains a non-analytic 
behavior as in the pure Ising case even for K4 = —0.5. 
This contrasts sharply with the lack of any traces of sin- 
gular behavior in the specific heat, cf. Fig. O From Fig. 
[T3I we see the same trend, namely that the susceptibility 
retains a non-analytic feature even for the largest K± val- 
ues we have considered, and which suffice to completely 
suppress the singularity in the specific heat. We have 
repeated these calculations with K xy — 0.3. The results 
are shown in Figs. [14] and [T5l with essentially the same 
results as in Figs. [T2l and [T3l 

We next attempt to estimate the critical exponents v 
and f3 for the model Eq. [1] for the set of parameters 
A = 1.0, K xy = 0.1, K A = -0.25. (For the same set of 
parameters, but K xy = and K4 = 0, see comments 
below). We base our calculations of these critical expo- 
nents on using the Binder cumulant Eq. [5] and the scaled 
staggered magnetization LP/ V M S , cf. Eq. [3 For these 
computations, we have used up to 3 • 10 6 sweeps over the 
lattice for each temperature. In addition, we have used 
Ferrenberg-Swendsen (FS) multihistogram reweighting 28 
of the raw data for the Binder cumulant in order to im- 
prove on the accuracy. The method of computation is de- 
scribed in Section Hi CI In Fig. [TH we show the Binder 
cumulant for various system sizes as a function of the 
temperature T without reweighting. The crossing points 
provide an estimate for T c . Even in the absence of FS 
reweighting, there is very little scatter in these cross- 
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FIG. 13: (Color online) The susceptibility of the staggered 
magnetization within each unit cell, Eq. [4] as a function of 
temperature T for the classical part of the generalized AT 
model Eq. [T] with A = 1.0 and K xy = 0.0, for various values 
of K 4 = 0.0,-0.1,-0.25,-0.5, and system size L = 512. 
Note that the susceptibility retains the non-analytical features 
of the Ising-case even for parameters where the specific heat 
anomaly is completely suppressed. 
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FIG. 14: (Color online) The staggered order parameter, Eq. 
|3l as a function of temperature T for the classical part of the 
generalized AT model Eq. [T] with A = 1.0 and K xy = 0.3, 
for various values of K4 = 0.0, —0.1, —0.25, —0.5, and system 
size L = 512. 
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FIG. 15: (Color online) The susceptibility of the staggered 
magnetization within each unit cell, Eq. [4] as a function of 
temperature T for the classical part of the generalized AT 
model Eq. [T] with A = 1.0 and K xy = 0.3, for various values 
of K 4 = 0.0, -0.1, -0.25, -0.5, and system size L = 512. Note 
the marked increase in the susceptibility as —K4 is increased, 
in contrast to the suppression of the anomaly in the specific 
heat. 
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FIG. 16: (Color online) The Binder-cumulant G, Eq. \5\ as 
a function of T for the model Eq. [T] for A = 1.0, K xy = 
0.1, K4, = —0.25 for various system sizes, in the absence of 
Ferrenberg-Swendsen reweighting. The inset shows a blowup 
of the temperature-region where the lines for various system 
sizes cross, providing an estimate for T c . 



ing points, and T c is determined with an uncertainty of 
much less than 1%. With Ferrenberg-Swendsen reweight- 
ing, this picture remains, as is seen from Fig. [T7] where 
reweighting is used. An accurate estimate for T c will turn 
out to be crucial in the following. Note also that the es- 
timates we get for T c from the crossing of lines in the 
Binder cumulant arc well in agreement from the some- 
what cruder estimates we would obtain from determin- 
ing the temperatures at which the peaks of the staggered 
susceptibilities occur. 

In Fig. 1181 we replot the same Binder-cumulant, now 
as a function of the quantity L X I V (T - T c )/T c , using es- 



timates for T c from Fig. [T7] and adjusting v to get data 
collapse. While we see that the above computations do 
not allow us to extract extremely precise values of v, it 
does allow us to conclude that the exponent v is con- 
sistent with the values obtained from the Ashkin- Teller 
model, and that v appears to be enhanced compared to 
the 2D Ising value v = 1. 

We next compute the quantity L@I V M S as a function of 
the quantity L 1 ^ {T — T c ) /T c to obtain the order param- 
eter exponent /3, by using the values of T c and v obtained 
from the scaled Binder cumulant in Fig. 1181 and then ad- 
justing (3 to get data collapse of all magnetization curves 
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FIG. 17: (Color online) The Binder-cumulant G, Eq. [5] as 
a function of T for the model Eq. □ for A = 1.0, if" = 
0.1, if4 = —0.25 for various system sizes, using Ferrenberg- 
Swendsen reweighting. The inset shows a blowup of the 
temperature-region where the lines for various system sizes 
cross, providing an estimate for T c . Note the consistency of 
the estimate for T c compared to Fig. 1161 
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FIG. 18: (Color online) The Binder-cumulant G, Eq. \5\ as a 
function of the quantity L 1/v (T - T c )/T c , for the model Eq. 
[[J for A = 1.0, K xv = 0.1, K 4 = -0.25 for various system 
sizes. Ferrenberg-Swendsen reweighting of the data is used. 
We have taken estimates for T c from Fig. [17] and adjusted 
the correlation length critical exponent v to achieve the best 
data collapse. As is seen, the optimal v is extremely sensitive 
to the chosen value of T c . 



for various values of L. The result of this procedure is 
shown in Fig. 1191 Again, from the above we cannot con- 
clude anything with great precision about the exponent 
(3, other than saying that it is consistent with the exact 
values that are known for the Ashkin- Teller model, i.e. 
Eq. [T] with K xy = 0. 

We have also checked the exponents for the same set 
of parameters as above, except that K xy = 0. We draw 
the conclusion that to the level of precision of the above 
computations, the exponents are not altered from the 
Ashkin- Teller case. However, when we repeat the pro- 
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FIG. 19: (Color online) The scaled staggered order param- 
eter L (itv M s , cf. Eq. O 

as a function of the quantity 
L 1/v (T - To)/T c , for the model Eq. □ for A = 1.0, K xy = 
0.1, iCt = —0.25 for various system sizes. Ferrenberg- 
Swendsen reweighting of the data is used. We have taken 
estimates for T c from Fig. [17] and estimates for v from Fig. 
1181 and adjusted the order-parameter exponent j3 to achieve 
the best data collapse. 

cedure for the same set of parameters as above, except 
that K xy = 0.3, we find that there is a clear deviation 
and that the exponents v and (3 definitely do not take 
Ashkin- Teller values. In particular, we get optimum data 
collapse for (3 clearly less then 1/8. From this, we infer 
that while the parameter K xy may be perturbatively ir- 
relevant, it may alter the universality class of the phase 
transition of the model if it is large enough. We also note 
that the reason that K xy appears to have much less of 
an effect on the transition when K4 = compared to 
when K4 — —0.25, is that the latter case represents a 
frustration of the ferromagnetic Ising ordering that low- 
ers the critical temperature of the system and enhances 
the effect of introducing K xy , which also frustrates the 
ferromagnetic Ising ordering, and promotes striped or- 
dering. 

Non-universality in (3 due to the presence of the pa- 
rameter K4 in the problem means that (3 in principle 
should vary slightly as we cross the pseudogap line verti- 
cally in the (x, T)-phase diagram of high-T c cuprates as 
the doping is varied, if we assume that the parameters 
of the effective model Eq. Q] varies as we move along the 
pseudogap line. In particular, a variation of [3 with K4 
is clearly seen from Fig. [T2J although we have not per- 
formed a detailed finite-size scaling analysis to determine 
(3 as a function of K4. We also note from Fig. 2] that in- 
troduction of K xy does not change the universality class 
of the transition when K4 = 0. We may therefore quite 
reasonably assume that the presence of K xy does not 
change the Ashkin- Teller universality class of the phase 
transition when K4 is present, provided K xy is not too 
large. We may then deduce that for negative K4, we will 
have -2 < a < 0, 1/8 < (3 < 1/4, and 7/4 < 7 < 7/2. A 
suppression of the specific heat anomaly as seen for the 
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case K4 = —0.25, puts us at a w —0.37, (3 w 0.15, and 
7 ps 2.07. The weak variation in the exponent from 
the Ising value 1/8 is due to the near-cancellation of the 
rather large, but opposite, variations in the specific-heat 
exponent a and the susceptibility-exponent 7, consistent 
with the scaling law a + 2(3 + 7 = 2. It is precisely 
the large variation in a that wipes out the specific-heat 
anomaly that also produces a large enhancement of the 
susceptibility of the staggered orbital magnetization, see 
Fig. El 



C. Comparison of Calculated Specific Heat with 
Experiments 



We use the results in Figs. [4lfTT1to estimate the peak 
value of the specific heat bump expected due to the tran- 
sition to see why it may be unobservable in experiments 
performed so far. In comparing with experiments, the fol- 
lowing should be borne in mind. The ordering below the 
transition temperature is three-dimensional. Hence, the 
observed specific heat will be that of the form calculated 
above as temperature is decreased towards T*, followed 
by a singularity characteristic of the 3D Ising model near 
T* and below it. However, the integrated value of specific 
heat divided by T under the singularity is only a fraction 
of the total entropy due to the loop order degrees of free- 
dom. As we discuss below, the latter itself is more than 
an order of magnitude smaller than the entropy due to 
fermionic excitations in the same temperature range. 

The area / (C v (T)/T)dT over all temperatures in each 
of the curves in Figs. I4TTT1 is 2 ln(2)/cB/unit-cell, reflect- 
ing that the calculations are performed for 2 Ising de- 
grees per unit-cell. Given that the ordered moment due 
to orbital currents is estimated in neutron scattering ex- 
periments to be 10 _1 /is/unit-cell, the integrated value is 
expected to be 2 ln(2)fcs/unit-cell multiplied by O(10 -2 ). 
To compare with experiments, we may consider calcula- 
tions for the case K xy = and \Ki/Ki\ between, say 0.25 
and 0.5. The peak value of the specific heat from Figs. 0J 
[TTlis then expected to be less than 0.5 x 10 _2 £;b per unit- 
cell or less than about 0.05 Joules/mole/degree. This 
should be compared with the measured specific heat 2 ^, 
which at about 200 K is about 200 Joules/mole/degree. 
In Ref. [29|, the electronic specific heat is estimated by 
subtracting the specific heat for a similar non-metallic 
compound to be about 2 Joules/mole/degree. Therefore, 
the bump has a peak which is 3 to 4 orders of magnitude 
smaller than the total specfic heat, and 1 to 2 orders of 
magnitude smaller than even the deduced electronic spe- 
cific heat. Given that the specific heat bump is spread 
out over temperatures of C(2T*), it is not surprising that 
with pseudogap temperatures of O(200) K or higher, it 
has gone undetected. There are underdoped cuprates 
with lower T* , in which a bump in the specific heat with 
magnitude of order that suggested here is claimed^ , to 
be observed. 



IV. UNIFORM SUSCEPTIBILITY 

Just as the onset of antifcrromagnetic spin-order has 
a weak parasitic non-analytic effect on the uniform mag- 
netic susceptibility, the onset of loop-current orbital mag- 
netic order may be expected to have a similar effect on 
the uniform magnetic susceptibility. Such an effect has 
indeed been measured recently in careful studies across 
T*(x)2.. 

Since the uniform magnetization is a parasitic effect on 
the staggered magnetization, we can calculate its tem- 
perature dependence by a Landau theory in which we 
consider the free energy for the staggered magnetization, 
but consider the minimal coupling of the uniform magne- 
tization to the staggered magnetization. Let M s be the 
staggered magnetization and (M) be the thermal aver- 
age of the uniform magnetization in the presence of an 
external field H . Let Fo(M s o) be the free-energy for the 
M s in the absence of an external magnetic field H. Quite 
generally, the leading terms in the free-energy are given 
by 



F (M S 



M 2 
2» 



- MH 



C 



-MzM 



(7) 



Here, C is a coefficient which gives the competition be- 
tween the staggered magnetization and uniform magne- 
tization. The sign of C is positive if as is reasonable, the 
staggered magnetization decreases if the uniform magne- 
tization increases, and vice versa. 

This form of the free energy gives correct answers only 
in the regime in which the staggered susceptibility is 
small and therefore is not valid very close to the tran- 
sition. Also, the susceptibility calculated is for magnetic 
field parallel to the direction of sub-lattice magnetization. 
In the simplest theory, this direction is perpendicular to 
the Cu-0 planes. In the experiments^, an angle closer 
to 7r/4 has been deduced for which some theoretical jus- 
tifications are provide d 16 ' 31 . Since the experiments are 
done in powder samples, we will ignore this issue for the 
present. 

Let Xso = (d 2 F /dM^ Q )~ 1 be the order parameter sus- 
ceptibility, which is calculated above. The subscript in 
Xso indicates the quantity in the absence of (M). xo is 
the uniform susceptibility in the absence of M s . Then in 
the presence of (M), induced by the external field H, the 
condition 



dF 
dM 



0. 



gives 



(M) 
Xo 



H + C{M){M*) = 



(8) 



(9) 



This gives (M) = xH in linear response (i.e. low H), 
with 



X 



Xo 



l + C X o(M 2 



(10) 
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Here, (M 2 ) is the thermodynamic squared magnetization 
in the presence of (M) , and x is the uniform susceptibil- 
ity. We may write Eq. [10] as 



X 



Xo 



(11) 



l + C* X o((M s > 2 +T Xs )' 
Also, quite generally, the order parameter susceptibility 



X s = 



d 2 F 

dM 2 



= X J + C<M 2 >= X J + CT X . (12) 



This gives 



Xs 



XsO 



1 + CTxXsO 



(13) 



Thus, using Eqs. (|12|13| . we may write x m terms of 
known quantities xo and Xos, to obtain 



CTxsoX + X - Xo = 0. 



(14) 



For T » T c , where the above treatment is valid, we 
have ACTxsoXo << 1 5 so that 



X « xo - CTxsoxl (T-T c )/T c »l. 



(15) 



The uniform susceptibility is therefore predicted to de- 
cline from its constant Pauli value at far above T c in 
the same range that XsO shows a rise. We suggest that 
the observed slow decrease of %(T) for temperatures well 
above T* be fitted to such a form. 

Well below T c , the model behaves as an Ising model. 
Therefore, the contribution of the ordered moments to 
the uniform susceptibility approaches zero exponentially 
as T -> 0. 



A. Mean-field Jump in dx/dT at T c 

In a mean-field calculation XsO does not change above 
T c . There is, however, a jump in dx/dT expected at T c . 
The experimental results have been quantified by such a 
jump 9 . To compare with available experimental results, 
we approximate Eq. (| 1 2[) as 



so that 



X ~ Xo - Cxi < M s > 2 , 



dx 2 d<M s > 2 

dT Xo dT 



(16) 



(17) 



Here, a temperature independent xo is assumed. We now 
need to know the right side of Eq. (fT7|) . This may be 
estimated as follows. Returning to Eq. ([7]), we may write 
F (M S ) as 

F (M S ) = &/2 {T T ^ ] M 2 + | Mt + .... (18) 



This defines the transition temperature T * in the absence 
of an external magnetic field H, i.e. for M = 0. It 
also defines an inverse susceptibility a for M s , which we 
expect to be of the same order as to the inverse of the 
density of states at the Fermi-surface, or equivalently of 
order Xo 1 - Combining Eq. ()18j) with the third term in 
Eq. |(7J), we see that a finite M leads to a decrease in the 
transition temperature 8T* , with 



ST* 



CM 2 /& 



Note also that 



M 2 «M S (0) 



2 (r * 



T) 



n 



(19) 



(20) 



where M s (0) 2 is the zero temperature value of M 2 . Us- 
ing this in Eq. (jXTJ) , the jump in the derivative of the 
susceptibility at T* is given by 



Hdx 
Xo dT 



CM s (0) 2 X o- 



(21) 



Now we need an estimate of CM s (0) 2 . This can be ob- 
tained from Eq. (|19p if we note that the transition tem- 
perature will be reduced to 0, i.e., = 1, for some 
magnetization M*. The magnitude of M* has to be the 
same order as M s (0) at zero field. Therefore 



C*M s (0) 2 /a« 1. 



(22) 



Using this above, the jump in dx/dT at T* is given by 



Tldx 
Xo dT 



«Xo 



(23) 



where we have used the estimate for a estimated earlier. 

In the experiments reported in Ref. H, a value of ^ 
between 0.2 and 0.3 has been deduced. This should be 
considered in good agreement with the estimate of 0(1). 
The weak assumptions in the analysis above are the lack 
of knowledge of the numerical constant between a and 
Xq^ 1 , and the unknown numerical constant on the right 
hand side of Eq. (|22|) . instead of 1. However, they cannot 
be off by more than an order of magnitude from those 
assumed. In a mean-field calculation %s0 does not change 
above T c . There is, however, a jump predicted in dx/dT 
expected at T c . 



V. SUMMARY 

We have studied the evolution of the specific heat and 
other thermodynamic properties in an effective theory 
of fluctuating orbital currents in high-T c cuprates. The 
motivation for the work has been to see if the finite- 
temperature break-up of a proposed ordering associated 
with a loop current pattern is consistent with both the 
existence of an order parameter in the pseudogap phase 
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below a temperature T*(x), and with an absence of an 
observed singularity in the specific heat and a weak singu- 
lar feature in the uniform magnetization at T*(x). This 
is a first step towards investigating, through quantum 
Monte Carlo simulations, whether the quantum break- 
up of such order gives rise to quantum critical fluctua- 
tions that could possibly explain the anomalous trans- 
port properties in the normal state of these compounds, 
as has been proposed in analytic calculations^. In this 
paper, we have shown that the effective field theory of 
the particular proposed order of orbital currents within a 
CuC>2 -plane passes this test by destroying the order while 
exhibiting no divergence in the specific heat. Instead, we 
have found bumps which we have estimated to be of a 
magnitude that are unobservable in experiments done so 



far. Moreover, we find a uniform magnetic susceptibility 
with a non-analytic behavior as a function of temperature 
as the phase transition is crossed. From a technical point 
of view, a principal result of our calculations is that the 
anisotropy considered in the Ashkin- Teller model as well 
as the next nearest neighbor interactions, in the range of 
parameters considered, are irrelevant perturbations. 
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